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Abstract 


A  three-dimensional  mathematical  model  of  a  solid  oxide  fuel  cell  is  presented,  which  allows  the  computation  of  the  local  distributions  of 
the  electrical  potential,  temperature  and  concentration  of  the  chemical  species.  The  physics  of  the  cell  and  the  simplifying  assumptions  are 
presented;  a  sketch  of  the  numerical  procedure  is  also  given.  The  numerical  results  obtained  with  hydrogen  as  the  fuel  are  compared  with 
results  from  other  simulation  codes  which  were  developed  for  a  planar  geometry.  The  numerical  results  show  the  behaviour  of  the  potential, 
temperature  and  current  distributions  when  certain  parameters  (geometry  of  the  cell,  electrolyte  materials,  temperature  in  the  channels)  are 
varied.  Numerical  simulation  is  also  used  to  obtain  an  optimum  for  some  geometry  parameters  such  as  cathode  thickness  or  rib  width. 
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1.  Introduction 

The  work  presented  here  deals  with  the  mathematical  mod¬ 
elling  and  numerical  simulation  of  natural  gas-fed  solid  oxide 
fuel  cells  (SOFCs)  under  a  stationary  regime.  The  principle 
of  an  SOFC  is  based  on  the  conversion  of  the  chemical 
energy,  which  is  stored  in  the  fuel  (hydrogen  or  methane), 
into  electrical  energy  through  an  electron-producing  electro¬ 
chemical  reaction  (see  e.g.  Ref.  [  1  ] ).  The  electrical  energy 
produced  by  an  electrochemical  reaction  depends  on  the  con¬ 
centration  of  the  reacting  species  through  the  Nemst  law.  A 
schematic  diagram  of  the  electrochemical  processes  is  given 
in  Fig.  1; 

1 .  oxygen  diffuses  through  the  porous  cathode  material; 

2.  oxygen  molecules  are  dissociated  and  ionized  at  the 
cathode/electrolyte  interface; 

3.  oxygen  ions  migrate  through  the  electrolyte  towards  the 
anode/electrolyte  interface; 

4.  fuel  diffuses  through  the  porous  anode  material;  and 

5.  hydrogen  contained  in  (and/or  produced  by)  the  fuel 
reacts  with  oxygen  ions,  producing  water  and  liberating 
electrons,  that  flow  back  to  the  cathode /electrolyte  inter¬ 
face,  via  an  external  circuit. 

Hydrogen  can  be  produced  from  methane,  and  an  interest¬ 
ing  feature  of  an  SOFC  is  the  fact  that  it  runs  at  a  temperature 
which  is  sufficiently  high  to  allow  internal  reforming,  which 
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Fig.  I .  Electrochemical  processes  in  a  SOFC. 


is  taken  into  account  in  the  present  model.  The  SOFC  systems 
seem  to  be  of  great  interest  for  the  use  of  natural  gas,  because 
of  their  high  power  generation  and  heat  recovery  efficiency, 
as  well  as  their  low  pollution  rate.  Moreover,  the  electrolyte 
in  an  SOFC  is  solid  rather  than  liquid  (as  for  instance  in  a 
phosphoric  acid  fuel  cell)  and,  therefore,  there  is  no  need  to 
handle  corrosive  liquids.  Because  of  their  high  temperature 
operating  conditions,  internal  reforming  of  methane  is  pos¬ 
sible,  along  with  heat  generation  from  the  hot  steam  which  is 
produced  by  the  cell.  However,  the  ideal  efficiency  rate  is 
never  attained  in  experimental  conditions.  Important  losses 
seem  to  be  internal  ohmic  losses;  diffusion  losses  and  over- 
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potential  due  to  the  electrochemical  reactions  can  also  be 
expected.  The  aim  of  the  mathematical  modelling  and  numer¬ 
ical  simulation  of  such  a  system  is  to  give,  for  a  given  set  of 
data  describing  the  geometry  of  the  cell  and  the  operating 
conditions,  an  estimate  of  these  losses  and  where  they  occur. 
Also,  by  varying  the  set  of  input  parameters  describing  either 
the  geometry  or  the  operating  conditions,  the  influence  of 
these  parameters  on  the  efficiency  of  the  cell  can  be  studied. 
Hence  modelling  may  be  viewed  as  a  design  tool. 

Numerical  simulation  of  SOFCs  has  already  been  per¬ 
formed  with  mathematical  models  which  were  obtained  with 
some  simplifying  assumptions  (see  Refs.  [2,3],  and  refer¬ 
ences  therein).  A  two-dimensional  mathematical  model  of 
the  cross  section  of  a  planar  SOFC  has  also  been  developed 
[4,5  ] ;  in  these  works,  as  in  the  present  work,  the  temperature, 
distribution,  concentration  of  the  species  and  electrical  cur¬ 
rent  are  obtained  within  the  unit  cell  itself,  whereas  in  other 
works,  which  are  concerned  with  the  modelling  of  a  whole 
stack,  the  solid  materials  were  assumed  to  be  all  at  the  same 
temperature,  and  the  current  density  was  computed  by  aver¬ 
aging  the  conductivities  of  the  various  components  [6-8] .  A 
number  of  models  of  the  unit  cell  also  exist  which  were  tested 
on  a  benchmark  defined  in  a  International  Energy  Agency 
( IEA )  workshop  [9-11].  The  originality  of  the  present  work 
lies  in  the  fact  that  it  yields  the  local  distribution  of  tem¬ 
perature,  electric  potential  and  concentration  in  a  three- 
dimensional  unit  cell  for  various  geometries;  it  uses  a  flux- 
conservative  discretization  scheme,  the  ‘finite  volume’ 
scheme  (seee.g.  Refs.  [12,13])  to  do  so. 

The  following  sections  are  organized  as  follows:  Section 
2  is  devoted  to  the  description  of  the  unit  cell  geometry  and, 
more  precisely,  to  the  geometries  which  have  already  been 
implemented  in  the  numerical  code.  Section  3  deals  with  the 
mathematical  model  itself,  i.e.  the  governing  equations  and 
the  boundary  conditions  which  were  defined  in  order  to  com¬ 
pute  the  local  distributions.  Section  4  describes  the  numerical 
procedure  and  the  computer  code  which  uses  this  procedure. 
The  validation  of  the  code  is  demonstrated  in  Section  5  by 
running  it  with  the  data  of  the  benchmark  defined  in  an  IEA 
workshop  [9]  and  comparing  it  with  other  numerical  codes. 
In  Section  6,  we  show  how  the  numerical  code  can  be  viewed 
as  a  design  tool.  First,  the  efficiencies  of  the  parallel  and 
cross-flow  cases  are  compared;  second,  the  influence  of  the 
electrode  thickness  on  the  efficiency  of  the  cell  is  studied.  In 
the  case  of  internal  reforming,  an  optimum  value  of  the  anode 
thickness  can  be  computed.  Finally,  numerical  results  for 
cylindrical  and  tubular  geometries  are  shown. 


2.  Description  of  the  geometry 

The  SOFC  system  is  in  fact  a  stack  of  electrochemical 
cells,  i.e.  it  is  made  up  of  several  repeating  ‘unit  cells’.  The 
layout  of  the  mathematical  equations  describing  the  physical 
and  chemical  phenomena  can  be  carried  out  on  the  smallest 
non-repeating  geometrical  pattern  which  represents  the  unit 


cell.  Several  concepts  exist  for  SOFCs,  based  on  the  geometry 
and  the  materials  used  for  the  different  components  of  the 
cell.  One  of  these  concepts  is  the  so-called  planar  bi-polar’ 
geometry,  see  Fig.  2.  This  concept  includes  the  ‘co-flow’  and 
‘counter-flow’  cases  where  the  air  and  gas  channels  are  par¬ 
allel,  and  the  ‘cross-flow’  case  where  the  channels  are  per¬ 
pendicular.  The  latter  case  is  preferred  by  industry  because 
of  the  simplicity  of  the  manifolding. 

The  results  of  Ref.  [3],  which  are  based  on  a  one-dimen¬ 
sional  model  along  the  channel  direction  (i.e.  a  two-dimen¬ 
sional  model  in  the  case  of  a  cross-flow  configuration) 
indicate  that  the  co-flow  geometry  should  be  the  most  effi¬ 
cient.  This  is  also  the  case  in  the  model  developed  in  Ref. 
[  1 4  ] .  In  a  previous  paper  [  1 5  ] ,  the  model  of  the  cross  section 
of  parallel  flow  configuration  was  developed,  ihis  model  is 
generalized  here  for  a  three-dimensional  unit  cell;  the  geom¬ 
etries  which  are  considered  are  the  planar  concept,  with  co-, 
counter-  or  cross-flow  designs;  our  numerical  code  can  also 
handle  the  tubular  concept  developed  by  Westinghouse, 
depicted  in  Fig.  3,  and  the  cylindrical  geometry,  sketched  in 


Fig.  2.  The  planar  concept  for  SOFCs. 
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Fig.  4.  The  code  is  also  designed  to  be  able  to  treat  other 
geometries  if  necessary.  The  unit-cell  mathematical  model  is 
independent  of  the  cell  geometry.  It  has  to  be  formulated 
before  carrying  out  the  numerical  computations  which  will 
yield  the  behaviour  of  the  physical  entities  for  one  unit  cell. 
It  can  also  be  used  to  compute  the  same  physical  entities  for 
a  whole  stack  provided  that  the  stack  is  made  up  of  the  same 
repeating  unit  cells,  and  that  the  boundary  conditions  at  the 
channel  walls  are  independent  of  the  cells  (i.e.  same  concen¬ 
tration  of  species  and  same  temperatures),  with  adequate 
boundary  conditions  on  the  exterior  of  the  stack  (thermal  and 
electrical  insulation).  So  as  to  leave  no  doubt,  the  various 
parameters  on  a  planar  bi-polar  co-  or  counter-flow  geometry 
are  given  in  Fig.  5.  The  cell  is  built  around  a  PEN  (positive 
electrode/electrolyte/negative  electrode)  sandwich,  where 
the  electrochemical  processes  occur.  The  interconnect  mate¬ 
rial,  I,  supports  the  structure  and  also  serves  as  the  current 
collector;  it  is  assumed  to  be  made  of  a  non-porous  material, 
and  the  diffusion  of  oxygen  only  takes  place  through  the 
cathode,  C.  The  fuel  methane  diffuses  through  the  porous 
anode,  A,  and  is  transformed  into  hydrogen  by  a  series  of 
chemical  reactions,  known  as  reforming,  in  the  presence  of  a 
catalyst,  throughout  the  anode.  A;  E  is  the  electrolyte,  and 
ns  the  solid  part  of  the  domain  of  study  (i.e.  anode,  cathode, 
electrolyte  and  interconnect);  Chf  and  Cha  are  the  fuel  and 
air  channels.  Note  that  in  a  two-dimensional  model  of  a  cross 
section,  the  temperature  and  concentration  are  assumed  to  be 
given  within  Chf  and  Cha,  and  therefore  the  domain  of  study 
is  fts  only  [  1 5  ] .  The  respective  boundaries  of  these  channels 
will  be  denoted  by  TChjl  and  rCh,- 
In  the  case  of  the  tubular  concept,  the  definition  of  the 
interfaces  must  be  somewhat  altered  to  take  into  account  the 
fact  that  the  electrolyte  is  insulated  from  the  interconnect;  the 
associated  interface  conditions  are,  however,  quite  easy  to 
implement. 


r» 

Fig.  5.  Geometry  of  a  unit  cell. 


3.  Mathematical  model  of  a  unit  cell 


3.1.  Conservation  laws 


Let  0  denote  the  electric  potential,  T  the  temperature  and 
X,  the  molar  concentration  of  species  i  ( i ' = 02,  N2,  H2,  H20, 
CO,  C02,  CH4).  The  potential  0  is  defined  throughout  fls, 
and  is  continuous  everywhere  except  at  the  electrode/elec¬ 
trolyte  interfaces  C/E  and  E/ A;  the  temperature,  T,  is  defined 
and  continuous  throughout  fts.  The  oxygen  molar  fraction 
X0,  (and  similarly  the  molar  fractions  X„  i=H2,  H20,  CO, 
C02,  CH4)  is  defined  and  continuous  in  the  cathode,  C  (or 
anode,  A).  The  mathematical  model  is  obtained  by  writing 
the  conservation  laws  for  the  heat  flux  q,  the  electric  current 
j,  and  the  mass  flux  Nh  i = 02,  N2.  H2,  H,0,  CO,  C02,  CH„: 

di v(/‘)  =0  in  I,  C,  E  and  A  ( 1 ) 

div(g)  =  Q  in  I,  C,  E,  A,  Cha  and  Chf  (2) 

div(fy)  =0  for  «= 02  and  N2  in  C  and  Ch.  (3) 

div(N,)=r, 

for  i=  H,,  H20,  CO,  CO,  and  CH4  in  A  and  Chf  (4) 

where  div  g = ijg,  /dx + 9g2/9y + dg3/dz  for  any  vector-valued 
function  g  =  (g„  g2,  g3)T,  Jt,  y,  z  representing  the  three  space 
coordinates.  Q  is  the  heat  source  term  (arising  from  the  pas¬ 
sage  of  electric  current  and  from  the  chemical  reactions),  and 
r,  the  production  rate  of  species  i.  The  fluxes  q,j,  (V,  and  the 
source  terms  Q  and  r,  must  now  be  related  to  the  unknowns 
T,  0  and  X,. 

In  the  above  equations,  it  is  assumed  that: 

1.  the  electron-chemical  reactions  occur  at  the  anode/ elec¬ 
trolyte  and  cathode/electrolyte  interfaces  (i.e.  heteroge¬ 
neous  reactions) ; 

2.  no  chemical  reactions  occur  within  the  air  channel,  and 

3.  reforming  and  water-shift  reactions  occur  within  the  fuel 
channel  and  anode  domains. 
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3.2.  Constitutive  laws 


both  the  porous  anode  and  the  fuel  channel.  The  reforming 
reaction: 


In  the  solid  parts,  the  thermal  flux  is  mainly  conductive 
and  modelled  by  Fourier’s  law  of  heat  conduction: 

q  =  grad  T  in  I,  A,  E  and  C  (5) 

where  grad  T  denotes  the  gradient  of  a  scalar  function  T,  i.e. 
grad  T=  (dT/dx,  dT/dy,  ar/3z)T,  and  A  is  the  thermal 
conductivity. 

In  the  gas  channels,  the  thermal  flux  is  mainly  convective 
in  the  gas  flow  direction,  and  conductive  from  the  channel  to 
the  solid  parts.  Hence  the  expression  for  the  flux  is: 
q=  £  Cpj  NjT—  A  grad  T  in  Chf  and  Cha  (6) 

where  G  is  the  set  of  various  components  of  the  gas  mixtures, 
i.e.  G=  {H2,  H20,  CO,  COa,  CH4}  in Chf,  whereas G=  {02, 
N2}  in  Cha.  The  electric  current  is  modelled  by  Ohm’s  law: 

j—-—<r  grad  0  in  I,  A,  E  and  C  (7) 

where  a  is  the  (temperature-dependent)  electric  conduc¬ 
tivity. 

In  the  porous  media  (i.e.  cathode  and  anode),  the  molar 
fluxes  are  mainly  due  to  diffusion  and  satisfy  Fick’s  law  of 
diffusion  in  a  mixture  ( see  e.g.  Ref.  [  16] ) : 

Ni  =  -  cD„„ uk  grad  X,+X,£ Nj  in  A  and  C  (8) 

where  c  is  the  concentration  and  Dumix  the  diffusion  coeffi¬ 
cient  of  species  i  in  the  considered  mixture.  In  the  gas  chan¬ 
nels,  the  molar  flux  is  mainly  convective  in  the  flow  direction 
and  due  to  molecular  diffusion  at  the  boundaries  rCh,  and 
rch.,  where  the  mass  exchange  may  be  taken  into  account 
using  a  mass-transfer  coefficient  (see  Sections  3.4  and  3.5). 
Hence  the  molar  flux  may  be  written  as: 

Nj  =  ~cVkX,  in  k  (for  *=Cha  and  Chf) 

for  i = 02,  N2,  H2,  H20,  CO,  C02  and  CH,  (9) 

where  Vk  is  the  velocity  of  the  mixture  in  channel  k.  In  fact, 
we  make  the  additional  assumption  that  the  nitrogen  flux  is 
zero  in  the  air  channel  and  porous  cathode  in  any  channel 
cross  section. 

3.3.  Source  terms 


CH4  +  H20  «■*  3H2 + CO  (11) 

is  irreversible  and  thus  the  heat  produced  is  equal  to  the 
enthalpy  A  HR.  Finally,  the  heat  source  term  can  be  expressed 


Q  =  cr  grad  0 grad  0+  Qchcm 

{Atfs  in  Chf 
AWs  +  A//Rin  A 
0  elsewhere 


(12) 

(13) 


The  volumetric  mass  source  term  r,  is  made  up  of  two 


n-aV+aV  (14) 

where  Is  and  r*  are  the  volumetric  reaction  rates  of  reactions 
(10)  and  (11),  respectively,  and  vs  is  the  stoichiometric 
coefficient  of  reactions  (10)  and  (11).  The  expression  for  r * 
was  taken  from  the  work  of  Lee  [  17] : 
i*  =  pk0  exp(  -  E^/RDP^P&ro  (15) 

where  p  is  the  resistivity  of  the  anode,  k0  the  frequency  factor, 
£R  the  activation  energy  the  reforming  reaction  and  a  and 
P  the  reaction  orders  of  CH4  and  H20.  The  expression  and 
the  data  given  by  Lee  [  17]  are  used  even  though  the  condi¬ 
tions  within  the  fuel  cell  anode  and  fuel  channel  are  different 
from  the  test  conditions.  The  H20  partial  pressure  level  and 
electric  field  which  are  present  in  the  SOFC  situation  [6] 
may  alter  the  reforming  rate.  Few  data  are  available  within 
the  open  literature. 

The  volumetric  reaction  rate  Is  is  computed  by  assuming 
the  shift  reaction  to  be  at  equilibrium: 


(16) 


y  XCOzXH; 

WHj0 

The  conversation  laws,  Eqs.(  l)-(4),  the  constitutive  laws, 
Eqs.  (5)-(9)  and  the  source  term  expression,  Eqs.  (12)- 
(14)  yield  the  following  system  of  (coupled)  elliptic 
equations: 

div(  -  <r  grad  0)  - 0  in  I,  C,  E  and  A  (17) 


div(  -  A  grad  T) 

=  a  grad  0  grad  0+  Qc ^  in  I,  C,  E  and  A  (18) 
div(  £  Cp,  N,T-  A  grad  7")  =  QctKm  in  Chf  and  Cha  ( 19) 


The  volumetric  heat  source  term,  Q,  is  in  fact  the  sum  of 
two  terms.  The  first  term  is  the  ohmic  heat  source  term  and 
can  be  written  tr  grad  0  grad  0.  The  second  term  is  the 
heat  produced  by  the  shift  and  reforming  reactions  which  take 
place  in  the  porous,  catalytic  anode,  and  we  shall  denote  it 
by  Gchem-  The  shift  reaction: 

C0  +  H20«C02  +  H2  (10) 

is  assumed  to  be  in  equilibrium.  The  heat  produced  is  equal 
to  the  term  A  Hs  of  the  reaction;  the  shift  reaction  occurs  in 


div(-cDo,mi,^^)=0inC  (20) 

div(  -cVch.X.)  =0  Tor  i=02  and  N2  in  CHa  (21) 

div(  -  cO„™,  grad  X.+X.X  Nj)  =  r, 

for  i= H2,  H20,  CO,  C02  and  CH,  in  A  (22) 

div(  —  cVchjX)  =r, 

for  i= H2,  H20,  CO,  C02  and  CH,  in  Chf 


(23) 
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3.4.  Interface  conditions 

Due  to  the  surface  electrochemical  reactions,  discontinui¬ 
ties  of  the  potential  and  heat  flux  at  the  electrode/electrolyte 
interfaces  exist.  Also,  in  the  three-dimensional  model,  the 
heat  and  mass  transfer  between  the  gas  channel  and  the  solid 
parts  will  be  taken  into  account  as  interface  conditions.  Note 
that  in  the  two-dimensional  model  these  become  boundary 
conditions,  since  the  temperature  and  the  concentration  in  the 
channels  are  then  data  to  the  model. 

The  interface  conditions  on  the  various  unknowns  and 
associated  fluxes  are  given  in  the  following  sections. 


#c-*E=^ln(X02)-/?o2[/-n]E  (29) 

-lF"i,n(^)"RH2[/',,lA  (30> 

where  AG  is  the  Gibbs  free  energy  change  of  reaction  (26), 
and  /?oj  and  J?Hj  are  the  polarization  resistivities. 

The  electric  current  is  continuous  everywhere  in  the  solid 
part,  and  is  linked  to  the  mass  flux  at  the  electrode/electrolyte 
interfaces  by  Faraday’s  law,  Eqs.  (27)  and  (28). 

At  the  channel/solid  interface,  the  electric  current  is  of 
course  zero,  hence  the  condition: 


3.4.1.  Temperature  and  heat  flux 

The  temperature  is  continuous  throughout  the  solid 
domain.  We  shall  model  the  boundary  layer  which  exists 
external  to  the  channel  walls  by  a  heat-transfer  law  of  the 
form: 

lqn]s  =  h(Ts-TG)  (24) 

where  [q  n]s  is  the  outward  normal  flux  to  the  solid  wall,  h 
the  heat-transfer  coefficient,  Ts  is  the  temperature  of  the  solid 
wall  and  Ta  the  temperature  of  the  gas. 

The  heat  flux  is  continuous  everywhere  except  at  the 
anode/electrolyte  interface,  where  the  electrochemical  reac¬ 
tion  generates  a  heat  source  term: 

[9-«]E+[*-n]A=-~^rASatE/A  (25) 

where  AS  is  the  entropy  of  the  reaction 
0.5O2  +  H2  «->  H20  (26) 

Note  that  this  model  does  not  take  into  account  the  effect 
of  radiation.  Work  is  under  progress  to  include  the  radiation 
equations  in  the  unit  cell  model  [18];  these  generate  a 
4th-order  nonlinearity  which  requires  a  specific  numerical 
procedure. 

3.4. 2.  Species  concentration  and  mass  flux 

The  species  concentrations,  X„  are  studied  in  the  channels 
and  the  porous  electrodes.  The  conservation  of  the  electrons 
must  be  written  at  the  interfaces  where  the  electrochemical 
reactions  occur.  Hence  the  mass  flux  is  related  to  the  electric 
current  by  Faraday’s  law: 

[(V0,-n]c=~jj~  at  C/E  (27) 

[AlHl-».]A=-^^atE/A  (28) 


3.4.3.  Potential  and  electric  currents 
The  potential  is  continuous  throughout  the  solid  parts 
except  at  the  electrode/electrolyte  interfaces  where  the  elec¬ 
trochemical  reactions  occur.  These  potential  jumps  are  mod¬ 
elled  by  Nemst’s  law  and  are  as  follows: 


<rgr2d4>n  =  0onrch,andrchf  (31) 

3.5.  External  boundary  conditions 

The  boundary  conditions  on  0  T  and  X,  (and/or  on  the 
associated  fluxes)  on  the  external  boundary  of  the  cell  remain 
to  be  stated.  These  boundary  conditions  depend  on  the  oper¬ 
ating  conditions  (see  Fig.  6).  The  boundary  conditions  con¬ 
sidered  here  are  based  on  the  idea  that  the  cell  is  part  of  a 
whole  stack  and  that  the  overall  behaviour  is  periodic.  There¬ 
fore,  on  the  lateral  boundary  IY  we  state  that: 

l/»]rL=0  (32) 

[?»lrt=0  (33) 

W-i «]rL=0  (34) 

for  «=02,  H2,  H20,  CO,  C02  and  CH„. 


TOTAL  FIXED  rOINT  M*TI AL  FIXED  KMNT 

Fig.  6.  Solution  procedure. 
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Periodicity  conditions  are  imposed  for  the  temperature  on 


the  top  and  bottom  of  the  cell,  whereas  a  me 
_/.lvg  is  imposed  for  the  current,  i.e.: 

:an  current  value 

m.-T=m,B 

(35) 

[9-n],-T=  -[9-nlt-B 

(36) 

[/•«].>=  [/«].„ 

(37) 

Lrn]rT=/.vg 

(38) 

Finally  the  temperature,  the  molar  fractions  and  the  molar 
flux  are  imposed  at  the  inlet  face  of  the  gas  channels. 


4.  Numerical  implementation 

4.1.  Numerical  procedure 

The  necessary  step  in  order  to  perform  a  numerical  simu¬ 
lation  is  to  discretize  the  set  of  equations  in  T,  <t>  and  X,.  The 
equations  themselves  are  of  the  Laplace  type  in  the  solid  part 
and  can  be  successfully  discretized  by  any  of  the  usual  meth¬ 
ods,  i.e.  finite  differences,  finite  elements  or  finite  volumes. 
However,  a  finite  volume  discretization  (see  Refs.  [  12,13] 
for  an  introduction)  was  preferred  over  the  finite  element 
method  because  it  is  computationally  cheaper  for  the  problem 
at  hand,  as  can  be  seen  from  the  results  shown  in  Ref.  [  15]. 
Moreover,  an  advantage  of  the  finite  volume  method  is  that 
it  is  ‘close  to  the  physics’  in  the  sense  that  it  is  based  on 
writing  the  balance  of  the  fluxes  through  the  boundaries  of  a 
discretization  cell  (which  is  called  a  ‘control  volume’),  and 
therefore  all  source  terms  can  be  accounted  for  in  a  natural 
way.  The  discretization  of  the  equations  yields  a  nonlinear 
finite  system,  where  the  unknown  parameters  are  the  average 
values  (over  the  control  volumes)  of  T,  <P and  Xh 

The  equations  in  T,  <P,  X,  are  coupled  by  the  electric  con¬ 
ductivity  (which  depends  on  the  temperature),  the  ohmic 
source  terms  and  the  interface  source  terms,  both  in  T  and  <P. 
A  rather  straightforward  iterative  method  is  used,  which  com¬ 


putes,  for  a  given  temperature,  the  solution  to  the  mass  and 
electric  problems  (X,  and  <P ),  and  once  convergence  is 
obtained,  solves  the  temperature  problem. 

4.2.  Computer  code 

The  computer  program,  called  SOFC3D,  consists  of  two 
parts.  The  first  part  is  the  computational  kernel  and  database 
manager.  It  was  written  using  the  C  language  and  contains 
1 1  000  instructions.  It  has  been  used  on  several  UNIX  sys¬ 
tems,  e.g.  Silicon  Graphics,  SUN  and  CRAY.  The  database 
input  file  contains  the  parameters  which  define  the  geometries 
(planar,  cylindrical  or  tubular),  the  composition  of  the  fuel 
(hydrogen  or  methane),  the  computational  mesh  and  the  set 
of  operating  conditions.  (Currently,  a  user-friendly  graphic 
X  interface  is  being  developed  [  19].)  The  second  part  is  a 
10  000  line  graphic  software  which  includes  a  user-friendly 
mesh  generator  and  a  three-dimensional  graphics  post¬ 
processor. 

The  results  which  are  presented  here  were  obtained  on  an 
INDIGO  SG  workstation  with  a  50  MHz  R4000  processor 
and  32  MB  of  real  memory.  The  database  file  requires  about 
100  MB  of  disk  memory. 

4.3.  Validation  of  the  numerical  code 

Since  we  have  no  experimental  data  to  compare  with,  we 
validate  our  code  by  comparing  the  results  obtained  for  a 
benchmark  which  was  defined  by  participants  of  the  IEA 
program  for  the  numerical  simulation  of  a  planar  geometry 
with  18  fuel  and  air  channels.  Note  that,  throughout  these 
tests,  in  accordance  with  the  IEA  benchmark,  the  porosity  is 
not  taken  into  account  and  the  thermal  conductivity  is  con¬ 
stant  and  equal  to  2  W  m~ 1  K“ '.  Here,  as  in  Ref.  [9],  the 
overpotential  at  each  interface  is  set  to  the  ohmic  loss  of  the 
electrolyte.  Results  are  given  in  Table  1  (results  from  IEA 
were  taken  from  Ref.  [9]).  It  appears  that  in  most  of  the 
models  cited  in  Table  1,  the  diffusion  effect  of  the  gaseous 
species  in  the  electrodes  is  only  accounted  for  in  the  direction 


Table  1 

Comparison  with  the  IEA  benchmark 


Max.  solid  temperature 


Min.  solid  temperature 


Domier(D)  0.689 

ECN  Petten  (NL) 

Eniriccrche  (I)  0.730 

Inst,  for Energiteknikk  (N)  0.71 

KFA-Julich  (D)  0.712 

Siemens  (D)  0.716 

Statoil  (N)  0.709 

Rise  (DK)  0.7101 

SOFC3D  ( ‘3D  diff )  0.7080 

SOFC3D  ( '  1 D  diff )  0.7203 


0.684  1085 

0.704 

0.722  1083 

0.7!  1084 

0.706  1073 

0.712  1062 

0.702  1082 

0.7034  1075 

0.6955  1055 

0.7079  1054 


1070  914  928 
1082  899 
1069  906  916 

1058  912  930 

1059  906  913 
1049  904  909 
1098  913  970 
1061  910  924 
1037  920  921 
1036  920  918 
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Table  2 

Gas  temperature  and  composition  at  inlet 


Gas  Temperature  (K)  02(%)  N2(%)  H2<%)  H20(%)  CO(%)  C02  <%)  CH,(%) 


Air  1173  21  79 

Fuel  I  1173  90  10 

Fuel  II  1173  26.26  49.34  2.94  4.36  17.10 


orthogonal  to  the  channel  axis.  We  implemented  both  this 
‘one-dimensional  diffusion'  and  the  three-dimensional  dif¬ 
fusion  effect.  Hence,  two  series  of  results  from  our  code  are 
listed  in  Table  1,  labelled  ‘ID  diff  and  ‘3D  diff .  The  ‘3D 
diff  cell  voltage  computed  with  our  code  is  somewhat  smaller 
than  the  average  cell  voltage.  However,  the  ‘ID  diff  cell 
voltage  which  we  compute  is  closer  to  the  voltage  obtained 
by  the  other  models.  We  show  in  Section  5  that  when  the 
three-dimensional  effect  is  accounted  for,  the  concentration 
profile  at  the  electrode /electrolyte  interface  is  not  constant 
and  that  the  low  values  of  concentrations  underneath  the  rib 
decreases  the  efficiency  of  the  cell.  The  effect  of  the  decrease 
in  activity  underneath  the  rib  is  recognized  to  be  also  of 
importance  in  stack  modelling  [20],  it  is  therefore  worth 
while  trying  to  measure  it  quantitatively  at  the  cell  level. 

The  code  can  simulate  a  cell  using  pure  hydrogen  only,  or 
using  methane  with  internal  reforming.  Tables  2  and  3  give 
the  general  operating  parameters  for  both  cases,  see  Ref. 
[21  ] .  These  parameters  will  be  used  throughout  the  numer¬ 
ical  runs,  unless  otherwise  mentioned.  The  following  section 
is  devoted  to  the  numerical  results  which  are  obtained  by 
running  the  code  on  different  geometries  of  adiabatic  cells, 
and  the  use  of  these  results  for  cell  optimization. 


5.  Using  the  code  as  a  design  tool 

5. 1.  Comparison  of  the  co-,  counter-  and  cross-flow  cases 

Numerical  tests  were  performed  in  the  case  of  the  co-  and 
counter-flow  for  a  planar  geometry  of  a  10  cm  long  unit  cell 
with  the  following  parameters:  active  area  5.42  mm  X 100 
mm,  anode  thickness  50  p,m,  cathode  thickness  50  pan,  elec¬ 
trolyte  thickness  150  pm,  channel  width  3  mm,  channel 
height  1  mm,  rib  width  2.42  mm,  and  total  bipolar  plate 
thickness  2.5  mm. 

Fig.  7  gives  a  comparison  of  the  efficiency,  i.e.  the  ratio  of 
electrical  energy  to  chemical  energy  at  inlet,  of  a  geometry 
with  5  fuel  and  air  channels  for  the  co-,  counter-  and  cross- 
flow  designs,  as  a  function  of  current  density.  In  this  numer¬ 
ical  test,  the  porosity  is  not  taken  into  account  and  the  thermal 
conductivity  is  constant  and  equal  to  2  W  m- 1  K_  1  as  in  the 
IEA  experiments.  The  overpotential  used  is  defined  in 
Ref.  [9]. 

These  computations  show  the  counter-flow  design  as  being 
the  most  efficient.  Note  that  this  is  in  agreement  with  the 
results  given  by  the  different  models  of  the  IEA  (see  Table 


1 )  but  it  contradicts  the  results  of  Ref.  [3],  where  only  the 
one-dimensional  effects  were  taken  into  account.  The  supe¬ 
riority  of  the  counter-flow  geometry  was  also  found  by 
numerical  experiments  on  a  three-dimensional  stack  [7] 
including  the  reforming  effect.  This  better  performance  can 
be  explained  by  the  fact  that  the  temperature  distribution  is 
better  for  the  efficiency  of  the  reforming  reaction. 

5.2.  Influence  of  rib  width  on  overall  efficiency 

The  conducting  material  between  two  channels  i: generally 
called  the  rib.  The  rib  width  seems  to  be  an  import;  . .t  param¬ 
eter  for  the  efficiency  of  the  cell: 

1 .  On  one  hand,  the  wider  the  rib,  the  smoother  is  the  elec¬ 
trical  current.  Hence  a  wider  rib  will  give  a  better  conduc¬ 
tion  of  the  electrical  current  and  less  ohmic  losses. 

2.  On  the  other  hand,  the  wider  the  rib,  the  n.  irrower  is  the 
channel;  when  the  channel  width  decreas;  >,  the  chemical 
species  do  not  diffuse  as  well  underneath  the  rib. 
Therefore,  the  numerical  code  can  be  used  to  check  if  an 

optimal  value  can  be  obtained  for  the  rib  width.  The  numerical 
results  confirm  the  competing  effects:  when  the  rib  width 
increases  (with  all  other  parameters  fixed),  the  ohmic  loss 
decreases  (Fig.  8),  but  the  potential  jump  also  decreases 
( Fig.  9).  With  all  other  parameters  set  to  the  values  indicated 
in  Tables  2  and  3,  an  optimal  rib  width  of  2.2  mm  is  obtained, 
as  indicated  in  Fig.  10. 

5.3.  Influence  of  electrode  thickness  on  efficiency 

We  study  here  the  influence  of  the  electrode  thickness  on 
the  concentration  of  H2  at  the  anode/electrolyte  interface  and 
on  the  concentration  of  02  at  the  cathode/electrolyte  inter¬ 
face.  It  is  clear  that  a  thicker  interface  electrode  will  be  more 
resistant  to  the  diffiision  of  the  gaseous  species;  in  order  to 
study  the  significance  of  this  effect  we  use  the  two-dimen¬ 
sional  pure  hydrogen  model  for  a  co-flow  geometry.  We 
found  that  the  cell  voltage  which  we  obtained  with  our  initial 
model  (labelled  ‘3D  difF)  could  be  significantly  increased 
if  the  diffusion  of  the  gas  in  the  porous  medium  was  taken  to 
be  only  in  the  direction  normal  to  the  electrolyte/anode  inter¬ 
face  (which  is  the  case  in  most  models  and  in  our  simplified 
model  labelled  ‘ID  diff).  In  our  model,  however,  the  dif¬ 
fusion  effect  is  accounted  for  in  all  directions,  and  the  hydro¬ 
gen  concentration  at  the  electrode/electrolyte  interface  is 
much  l^v'er  away  from  the  channel  than  close  to  it.  This  is 
clear;.,  seen  in  Fig.  11,  which  represents  the  hydrogen  con- 
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Table  3 

General  operating  conditions 

Average  current  density 

Fuel  utilization 

Overpotential 

300  mA  cm-2 

85% 

7 

see  Ref.  [8] 

Electrical  <r  (fi1  rrT  ')  and  thermal  A  ( W  m 

1  K- ')  conductivities 

Cathode 

Electrolyte 

Anode 

Interconnect 

a  4.2  X107 

Jim)  3  34X  10*exi/— — — — 1 

9.5  X107  J  1150) 

9.3  X  10s  J  1100) 

T 

A  50% 

A  3 

1^  r  |  J.J4X  lu  exp^  ^  j 

r- expl  Tj 

50% 

3 

T  6XPl  T  ) 

3.5 

Diffusion  coefficient 

Dn 

Oh, 

D,  (Fuel  II) 

Values  (cm2  s'1) 

0753X®'5 

omx\£j‘ 

Fig.  7.  Efficiency  of  co-,  counter-  and  cross-flow  geometries. 

centration  profile  for  various  anode  thicknesses,  using  pure 
hydrogen  as  a  fuel. 

Fig.  12  represents  the  PEN  potential  jump  computed  with 
the  Nemst  law  with  respect  to  the  molar  fractions  in  the 
channels,  at  the  electrode/electrolyte  interface,  and  with  the 
maximum  value  of  the  molar  ff action  at  the  electrode  /elec¬ 
trolyte  interface.  Taking  this  last  value  is  in  fact  assuming  a 
one-dimensional  diffusion  in  the  porous  electrode.  It  is  clear 
from  Fig.  1 2  that  the  three-dimensional  effect  in  the  diffusion 
loss  is  quite  important,  and  that  it  can  in  fact  become  dramatic 
in  the  case  of  a  narrow  anode. 

In  the  case  where  methane  is  used,  the  above  results  will 
be  altered  since  the  methane  reforming  rate  increases  with 


the  anode  thickness  (and  the  amount  of  catalyst),  assuming 
a  volumetric  reforming  reaction.  There  are  now  several  com¬ 
peting  effects  which  play  a  role  in  the  optimization  of  the 
anode  thickness: 

1.  The  diffusion  of  gas  is  better  for  a  relatively  thin  anode 
(but  not  too  thin,  as  observed  in  Fig.  12) . 
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800  pm. 


Fig.  12.  PEN  potential  jump  vs.  anode  thickness. 

2.  Tits  production  of  methane  by  volumetric  reforming 
increases  with  the  thickness  of  the  anode  (see  Fig.  13). 

3.  The  methane  reforming  is  an  endothermic  reaction,  so  that 
when  the  methane  production  increases,  it  decreases  the 
temperature. 

4.  The  ohmic  loss  increases  with  the  anode  thickness,  which 
tends  to  increase  the  temperature. 

The  numerical  code  simulating  the  three-dimensional 
model  for  methane  in  the  planar  geometry  was  used  to  study 
this  phenomenon.  The  input  parameters  were  set  as  in  Tables 
2  and  3,  with  a  fixed  inlet  fuel  flux  (0.038“'  mol/h).  The 
thickness  of  the  anode  was  varied  from  50  to  800  pm.  In  fact 
for  an  anode  thickness  of  50  pm,  there  is  not  sufficient  hydro- 


Fig.  IS.  Nemst  potential  at  the  electrolyte  level  vs.  anode  thickness. 


gen  produced  by  reforming  for  the  requested  produced  cur¬ 
rent,  and  the  numerical  code  cannot  reach  a  solution. 

The  results  (Fig.  14)  show  that,  for  a  planar  geometry,  a 
thickness  of  200  pm  seems  to  be  optimal  with  respect  to  the 
efficiency  ( defined  as  the  ratio  of  the  electric  power  over  the 
chemical  power  of  the  inlet  gas).  It  can  be  observed  from 
Fig.  15  that  the  Nemst  potential  increases  steadily  with  the 
anode  thickness,  but  so  do  the  ohmic  losses  (Fig.  16)  and  the 
polarization  losses  (Fig.  17).  Note  that  the  fuel  utilization 
decreases  with  the  anode  thickness  (Fig.  18). 
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Fig.  17.  Polarization  losses  vs.  anode  thickness. 


The  minimal  temperature  may  go  below  the  inlet  temper¬ 
ature  of  the  fuel  because  of  the  endothermic  characteristic  of 
the  reforming  reaction,  as  may  be  seen  in  Fig.  19.  This  is  in 
agreement  with  the  remarks  on  a  ‘cold  cell’  of  Ref.  [22],  The 
rib  effect  may  again  be  seen  in  Figs.  20-22. 

5.4.  Results  for  other  geometries 

The  adaptivity  of  our  three-dimensional  SOFC  code  can 
be  demonstrated  by  running  numerical  simulations  on  other 
fundamental  geometries.  A  sample  of  results  concerning  a 
tubular  (Westinghouse  concept)  geometry  and  a  cylindrical 
geometry  are  given. 
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5.5.  Tubular  geometry 

In  this  test,  the  mean  current  density  is  taken  equal  to  1500 
mA  cm'2  in  the  interconnect.  The  dimensions  of  the  cell 
were  taken  from  Refs.  [  23,24  ]  and  are  shown  in  Fig.  23.  This 
figure  also  shows  the  contours  of  current  density  at  the  outlet 
section  of  the  cell.  It  can  be  noted  that  the  highest  density  is 
near  the  intcrconnect/anode  interface;  this  is  due  to  the  fact 
that  the  anode  is  thin,  and  that  the  interconnect  area  is  not 


very  large.  On  the  cathode  side,  ho  wever,  the  current  density 
is  lower  because  the  cathode  is  thicker.  Fig.  23  shows  the 
contours  of  the  molar  fraction  of  oxygen  and  the  temperature 
at  the  outlet  section  of  the  cells.  Note  that  the  minimum  of 
the  molar  fraction  of  02  and  the  maximum  temperature  are 
reached  in  the  region  of  highest  current  density  (see  Fig.  24). 
Fig.  25  shows  the  comparison  of  the  potential  drops  for  a 
planar  geometry  and  a  tubular  geometry  which  were  chosen 
such  that  the  input  electric  current  is  the  same  on  the  upper 
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Fig.  25.  Comparison  of  the  tubular  and  planar  geometries. 


surface  of  both  geometries.  These  results  show  that  the  tubu¬ 
lar  geometry  is  more  efficient. 

5.6.  Cylindrical  geometry 

In  the  absence  of  data  for  a  so-called  ‘cylindrical’  geom¬ 
etry,  we  derived  the  dimensions  of  a  cylindrical  cell  by  dis¬ 
torting  a  planar  cell.  The  dimensions  of  the  cell  are  shown  in 
Fig.  26.  The  computations  were  performed  on  a  ‘slice*  of  the 
cell  only,  taking  symmetry  considerations  into  account.  The 
mean  current  density  is  taken  equal  to  1 80  m  A  cm  "  2.  Fig.  27 
shows  the  distribution  of  temperature  in  the  electrolyte  and 
Fig.  28  the  distribution  of  hydrogen  in  the  anode.  Again,  one 


can  see  that  the  diffusion  is  always  weaker  at  the  centre  of 
the  rib. 


6.  Conclusions 

A  three-dimensional  numerical  code  was  developed  for  the 
numerical  simulation  of  SOFCs  and  stacks,  which  can  be 
used  for  various  geometries  of  SOFCs.  The  mathematical 
model  is  based  on  a  local  balance  approach,  and  the  finite 
volume  method  was  found  to  be  well  adapted  for  our  problem. 
N  umerical  runs  were  performed  first  to  validate  our  code  with 
a  benchmark  which  was  defined  by  the  participants  of  an  IGA 
workshop;  the  results  which  were  obtained  with  our  code 
were  in  good  agreement  with  those  of  the  IEA  workshop  for 
the  test  case  (with  hydrogen  as  a  fuel),  except  for  the  fact 
that  our  code  seemed  to  be  a  little  more  pessimistic  than  most 
others  with  respect  to  the  electrical  potential  of  the  cell.  This 
is  probably  due  to  the  diffusion  effect  which  our  model  takes 
into  account  more  accurately.  The  code  was  then  used  as  a 
design  tool  for  the  cell  geometry;  we  compared  the  cross-, 
counter-  and  co-flow  designs  for  a  planar  geometry  and  found 
that  the  counter-flow  was  optimal  (with  hydrogen  as  a  fuel) . 
We  also  found  that  we  could  use  the  code  to  compute  the 
optimal  value  of  one  given  parameter,  such  as  the  rib  width, 
or  the  anooe  thickness.  In  the  latter  case,  the  reforming  effect 
was  studied  (using  methane  as  a  fuel ) .  Comparisons  between 
the  tubular  and  the  planar  geometry  were  carried  out,  and 
showed  less  ohmic  loss  for  the  first  geometry.  The  flexibility 
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of  our  code  allows  the  simulation  of  several  geometries. 
Indeed,  a  new  geometry  could  be  introduced  relatively  easily 
thanks  to  the  modularity  of  the  code. 

Further  work  includes  the  study  of  the  effect  of  the  radia¬ 
tion  inside  the  fuel  and  air  channels.  In  Ref.  [  22] ,  the  effect 
of  radiation  was  introduced  at  the  stack  level  by  means  of  the 
computation  of  an  ‘effective  thermal  conductivity’  [22] ,  and 
it  was  shown  in  Ref.  [25],  using  a  finite  element  approach 
for  temperature  and  electrical  current,  that  the  radiation  effect 
has  a  definite  influence  on  the  temperature  distribution.  A 
sensitivity  analysis  with  respect  to  the  geometry  parameters, 
the  physical  properties  of  the  constitutive  materials  and  the 
operating  conditions  is  also  planned. 

7.  List  of  symbols 
A  anode 

c  concentration,  mol  m  ~ 3 

C  cathode 
Cha  air  channel 

Chf  fuel  channel 

Cpi  heat  capacity  of  gas  i,  kJ  kmol  '  K“ 1 
A,™*  diffusion  coefficient  of  species  i  in  mixture 
E  electrolyte 

£  activation  energy,  kJ  kg  1  K " 1 
F  Faraday  constant,  A  s  mol  ” 1 
h  heat  transfer  coefficient,  W  m  2  K"  1 
j  electric  current,  mA  cm  ~  2 

K  equilibrium  constant 

Nj  molar  flux  for  species  i,  mol  s  “ 1 
q  heat  flux,  Ws“‘ 

T  temperature,  K 

Q  heat  source  term,  W  m  “  3 

Cchem  heat  source  term  due  to  chemical  reaction,  W  m  _  3 

V  velocity  of  gas  in  channels,  ms'1 

Xj  molar  concentration  of  species  i 

Greek  symbols 

AW  enthalpy,  kJ 

AS  entropy,  kJ  K~ 1 

Teh,  interface  between  air  channel  and  solid  walls 
rChf  interface  between  fuel  channel  and  solid  walls 
A  heat  conductivity,  W  m“ 1  K~ 1 

0  electrical  potential,  V 

p  electrical  resistivity,  W  ~ 1  m  K 

a  electrical  conductivity,  W  m  " 1  K  ” 1 

fts  solid  part  of  the  physical  domain 

Subscripts 


f  fuel 

i  chemical  species 

R  reforming  reaction 

S  shift  reaction 
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